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Abstract 


Understanding  the  concepts  of  spectral  estimation  may  become  increasingly  important  in  future 
radar  system  design.  The  purpose  of  this  technical  memo  is  to  introduce  these  concepts.  A  simple 
model  is  formulated  using  a  uniformly  spaced  linear  phased  array,  and  the  classical  Fourier, 
Maximum  Entropy  (ME),  and  Multiple  Signal  Classification  (MUSIC)  techniques  are  described. 
The  performance  of  these  three  techniques  is  discussed  through  qualitative  comparison. 
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T.  Introduction 


Many  spectral  estimation  methods  have  been  developed  to  determine  the  spatial  distribution 
of  sources  in  the  object  field  of  an  antenna  system.  These  techniques  are  useful  for  both  resolving 
the  presence  of  closely  spaced  sources,  and  detecting  individual  point  scatterers  on  a  single  source 
for  imaging  applications.  The  intent  of  this  technical  memo  is  to  provide  the  necessary  background 
to  understand  these  techniques.  The  presentation  begins  with  the  development  of  a  simple  model 
using  a  uniform  linear  phased  array  as  a  basis  for  the  concepts.  Three  different  spectral  estimation 
techniques  are  discussed:  the  classical  Fourier  method,  the  maximum  entropy(ME)  method,  and  the 
multiple  signal  classification(MUSIC)  method.  The  Fourier  method  is  often  termed  conventional 
and  non-parametric,  while  the  ME  and  MUSIC  methods  are  parametric  and  typically  referred  to  as 
superresolution  techniques.  The  performances  of  these  techniques  are  compared  qualitatively  with 
respect  to  their  underlying  concepts,  with  reference  to  the  open  literature.  The  presentation  is 
restricted  to  the  narrowband  case  with  all  signals  having  the  same  earner  frequency. 


n.  Model  Formulation 


Consider  the  linear  phased  array  shown  in  figure  1.  The  arriving  wavefront  at  the  aperture 
is  the  resultant  sum  of  M  radiating  sources,  and  is  sampled  by  each  of  the  N  elements.  The 
complex  amplitude  of  the  noiseless  signal  s(n)  at  each  element  may  be  expressed  as 

M 

(1)  s(n)  =  exp[jkd(n-^)sin((p j] 

m=l 

where  N  is  the  total  number  of  elements,  M  is  the  total  number  of  arriving  signals,  k  =  2%fk  is  the 
wave  number,  d  is  the  element  spacing,  Am  is  the  complex  amplitude  of  the  m*  signal,  and  is  the 
angle  of  arrival  of  the  m*  signal.  The  array  element  locations  are  referenced  to  the  center  of  the 
array  as  depicted  in  figure  1 . 
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figure  1. 

The  signal  measured  at  each  array  element  is  perturbed  by  additive  noise.  The  noise  is  principally 
introduced  by  the  antenna  element  and  associated  receive  module,  and  is  normally  characterized  by 
their  combined  noise  figme.  In  practice,  the  noise  is  modeled  as  a  white,  stationary,  ergodic 
random  process  [2],[3].  For  the  narrow  band  case,  the  complex  amplitude  of  the  noise  signal 
present  at  each  element  is 

(2)  w(n)  =  expHpJ 

where  the  amplitude  Bn  is  Rayleigh  distributed,  and  the  phase  Pn  is  uniformly  distributed  over 
(0,27c).  The  complex  signal  w(n)  is  Gaussian  distributed,  normally  with  zero  mean.  Combining 
equations  (1)  and  (2),  the  observed  signal  at  each  element  is 

(3)  x(n)  =  s(n)  +  w(n) 

The  following  sections  are  based  on  the  simple  model  presented  in  this  section. 


ITT.  Classical  Fourier  Approach  to  Spectral  Estimation 


Consider  the  noise  free  expression  for  the  signal  at  each  element  given  in  equation  (1).  An 
interesting  base-line  approach  may  be  taken  to  estimate  the  relative  strength  and  angle  of  arrival  of 


the  in^  signal.  Recalling  a  common  approach  used  in  antenna  system  design,  the  far  field 
expression  of  the  aperture  radiation  pattern  may  be  determined  by  the  Fourier  transform  of  the 
current  distribution  across  the  aperture.  In  this  case,  the  far  field  source  distribution  may  be 
determined  from  the  inverse  Fourier  transform  of  the  aperture  excitation  current.  That  is. 


(4) 


L/2 

f(u)  =  I  i(x)  exp[-jkxu]  dx 
•'-L/2 


where  f(u)  is  the  normalized  complex  far  field  source  distribution,  i(x)  is  the  normalized  complex 
aperture  excitation  current,  L  is  the  total  length  of  the  aperture,  and  u  =  sin((p).  For  a  uniformly 
spaced  linear  array,  the  aperture  excitation  may  be  expressed  as  a  uniformly  sampled  continuous 
aperture  excitation.  Thus, 


(5) 


N 


i(x)  =  i^(x)  ^  5(^x-d(^n-^^]) 

n=l 


where  ic(x)  is  the  continuous  current  excitation.  Combining  equations  (4)  and  (5)  results  in 

N 

(6)  f(u)  =  ^  i(n)  expj^-jkd(n-^^)u 

n=l 

Observing  equations  (1)  and  (6),  it  is  apparent  that  s(n)  is  proportional  to  i(n).  Substituting  s(n) 
into  equation  (6),  and  considering  only  the  m*  signal  gives 

N 

(7)  f(u)  =  ^  exp[-jkd(n-^)(u-u^) 

n=l 

where  u  =  sinfffi  ).  Note  that  equation  (7)  is  a  geometric  series,  and  may  be  simpUfied  to 

sin 

(8)  f(u)  =  A^  — 


N 


in 


u-u 


m 


=  A  f  (u-u  ) 
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where  f<,(u-Un,)  represents  the  radiation  pattern  of  a  uniformly  illuminated  linear  array,  whose  beam 
may  be  steered  to  a  position  9^  by  introducing  a  linear  phase  shift  across  the  aperture.  A  plot  of 
fo(u-Unj)  with  Unj=0,  is  shown  in  figure  2. 


f(u) 


The  first  sidelobe  is  -13dB  relative  to  the  main  lobe,  and  successive  sidelobes  fall  off  at  6dB  per 
octive.  This  pattern  illustrates  many  of  the  problems  associated  with  spatial  processing,  and  in 
particular  those  pertaining  to  the  performance  of  the  Fourier  technique.  It  is  apparent  that  a  main 
lobe  will  occur  in  the  direction  of  a  defined  angle  of  arrival,  but  ambiguities  exist  in  the  directions 
of  the  sidelobes.  In  addition,  the  resolving  power  is  determined  by  the  width  of  the  main  lobe. 
When  two  or  more  sources  are  present,  superposition  may  be  applied  to  express  the  angle 
spectrum  of  the  total  field  present  at  the  aperture. 


M 


(9) 


F(u)  =  X 


sin  N — ^ 
V  2 


m  /  u-U 
m=l  .  I  m 


Sin 


Equation  (9)  represents  M  complex  weighted  radiation  patterns  positioned  throughout  the  angle 
spectrum  corresponding  to  M  discrete  angles  of  arrival.  Consider  the  simple  case  of  two  returns 
whose  radiation  patterns  overlap.  It  is  possible  that  the  presence  of  both  returns  may  not  be 
detected,  depending  upon  their  relative  strength  and  phase  alignment.  Thus,  in  the  case  of 
resolving  multiple  returns  having  overlapping  radiation  patterns,  the  Fourier  approach  appears  to 
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have  considerable  limitations.  One  way  to  improve  the  performance  is  to  taper  the  aperture 
distribution^].  This  may  reduce  the  sidelobe  level  considerably,  at  the  expense  of  widening  the 
main  lobe.  There  are  a  variety  of  aperture  tapers  presented  in  the  literature  that  realize  this 
exchange[3].  The  effective  performance  may  be  extended  further  by  adaptively  controlling  the 
aperture  taper.  Thus,  the  trade-off  between  the  sidelobe  level  and  the  width  of  the  main  lobe  could 
adapt  to  the  situation.  The  resolving  power  of  the  classical  Fourier  approach  is  at  best  on  the  order 
of  a  standard  3dB  beamwidth(VL),  depending  upon  the  relative  source  strength  and  phase 
alignment,  and  the  SNR[3]. 

Despite  the  limitations  presented,  there  are  some  notable  advantages  to  the  Fourier 
approach.  First,  it  is  non-parametric  in  the  sense  that  no  known  information  of  the  received  signal 
is  necessary  for  implementation.  Second,  it  is  robust  because  it  is  relatively  insensitive  to 
parameter  changes.  And  third,  it  is  easy  to  implement  making  it  a  practical  approach. 


TV.  Parametric  Spectral  Estimation 

Consider  the  principles  behind  the  basic  Fourier  approach  just  described.  It  is  clear  that  its 
limitations  result  from  the  assumption  that  the  field  is  zero  beyond  the  aperture.  The  principal  idea 
behind  parametric  spectral  estimation  (superresolution)  is  to  not  make  this  assumption. 
Alternatively,  the  signal  information  available  at  the  aperture  is  used  to  estimate  additional  signal 
information  beyond  the  aperture.  The  pioneering  work  in  this  area  was  performed  by  Burg  [5]  in 
1967,  and  is  commonly  known  as  the  Maximum  Entropy  Method(MEM).  The  maximum 
entropy(ME)  approach  appears  throughout  the  literature  as  the  basis  for  solving  many  spectral 
estimation  problems.  Since  its  inception,  a  number  of  superresolution  techniques  have  been 
developed.  Perhaps  the  most  popular  being  the  Multiple  Signal  Classification  (MUSIC)  algorithm 
developed  by  Schmidt[9]  in  1979.  As  with  most  parametric  spectral  estimation  techniques,  both 
MUSIC  and  ME  operate  on  the  spatial  correlation  or  covariance  matrix  of  the  observed  signal.  The 
covariance  matrix  may  be  developed  by  first  summing  equation  (3)  over  N,  and  expressing  it  in  a 
more  compact  form  using  matrix  notation.  That  is, 

(10)  x=s+w  =  Ba-t-w 
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The  observed  signal  vector  x  represents  a  single  snapshot  of  data  corresponding  to  a  particular 
instant  of  time.  In  general,  data  is  collected  over  a  short  period  of  time  such  that  x  consists  of  a 
number  of  independent  snapshots.  Therefore  the  measured  data  normally  available  for  processing 
is 
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(16) 


x(p)  =  Ba(p)  +  w(p)  p- 1,2,...P 


where  P  is  the  total  number  of  snapshots.  The  angles  of  arrival  are  assumed  not  to  change 
significantly  in  the  processing  interval,  which  is  shown  in  equation  (16)  where  the  direction  matrix 
B  is  not  a  function  of  p.  The  signal  in  space  vector  S  however,  is  a  function  of  p  because  sources 
responsible  for  its  generation  are  generally  unpredictable,  and  therefore,  modeled  as  stochastic 
processes!  1].  Under  this  assumption,  temporal  averaging  may  be  performed  to  improve  the 
spectral  estimation. 

The  spatial  correlation  of  the  observed  signal  vector  x  is  defined  as 

(17)  R  =  E[x(p)x'^(p)] 

where  E[-]  is  the  expectation  operator,  and  +  denotes  complex  conjugate  transpose.  Substituting 
equation  (16)  into  (17)  and  noting  that  the  noise  and  signal  are  statistically  independent  results  in 

(18)  R  =  E[Baa'^B'^]  +  E[ww'^] 

=  brX+r,=  r,+r, 

where  Ra,  R^  and  Rs  are  the  spatial  correlation  matrices  of  the  signal  in  space  vector,  noise 
vector,  and  received  signal  vector,  respectively.  If  the  noise  is  white  with  zero  mean  and  variance 
o^,  the  spatial  correlation  matrix  of  the  noise  vector  may  be  expressed  as 

(19) 

and  equation  (18)  becomes 

(20)  R  =  Rs 

where  R  is  a  square  matrix  whose  rank  corresponds  to  the  number  of  elements,  and  the  total 
number  of  sources  that  may  be  detected  is  N-1.  The  mathematics  are  also  valid  for  the  single 
snapshot  case. 
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A.  Maximum  Entropy  (ME)  Method 


The  principle  behind  the  ME  approach  is  to  choose  the  spectmm  f(u)  whose  inverse  Fourier 
transform  is  i(x)  for  Ixl  <  L/2,  while  making  the  fewest  assumptions  regarding  unmeasured 
information  in  the  region  1x1  >  L/2.  This  condition  is  the  information  theory  concept  for 
maximizing  entropy.  From  the  work  of  Shannon  in  1948,  entropy  was  found  to  be  proportional  to 
the  logarithm  of  the  power  spectrum,  and  is  commonly  expressed  as 

oo 

(21)  H  =  j  ln[S(u)]du 


where  S(u)  =  IF(u)|2.  An  extension  of  the  classical  Fourier  approach  described  in  the  previous 
section  may  be  thought  of  in  terms  of  the  Wiener-Khintchine  theorem.  This  theorem  states  that  the 
power  spectral  density  and  the  autocorrelation  function  are  Fourier  transform  pairs.  Using  the 
notation  of  equation  (21),  this  may  be  stated  as 


(22)  R(±x^)  =  J  S(u)  exp  [±  jkx^u]  du  n=0,l,...,(N— 1) 


where  R(Xjj)  is  the  autocorrelation  function  of  the  field  across  the  array.  Mathematically,  the  ME 
method  consists  of  maximizing  equation  (21),  while  satisfying  equation  (22).  For  the  case  of  the 
periodic  array,  the  ME  estimate  of  the  power  spectrum  may  be  expressed  in  closed  form  [4]  as 


(23) 


'MEM 


(U)^ 


Tl 


|[R"Ww(u)P 


where  dj^  =  (1,1)  element  of  R'^,  5  =  [1,0,0,...,0]’^,  and  W(u)  is  the  direction  of  look  vector 
referenced  from  the  center  of  the  array.  That  is. 


(24) 


W(u)  = 


exp[jkdu(l-N/2] 


exp[jkdui 


3-n)/2] 


exp[jkdu(N-l)/2] 
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The  resolution  performance  of  the  ME  approach  when  the  true  autocorrelation  function  is 
known  is  noise  limited,  and  is  supperior  to  the  Fourier  approach  for  both  coherent  and  incoherent 
sources.  It  is  shown  in  reference  [8]  that  the  resolution  limit  expressed  in  fractional  beamwidths  is 
pj-oportional  to  where  0. 1X,/L  resolution  requires  ~  33dB  of  SNR.  This  hmitation 

has  been  found  to  be  a  result  of  uncompensated  noise  bias  present  in  the  diagonal  elements  of  the 
autocorrelation  matrix[6]  (see  equations  (18)  thru  (20)).  If  the  effect  of  the  noise  bias  is 
compensated,  angular  resolution  approaching  0.1  X/L  may  be  achieved  at  much  lower  values  of 
SNR[6],[7].  Thus,  the  ME  approach,  with  or  without  noise  bias  compensation,  will  out-perform 
the  classical  Fourier  approach  for  the  case  of  incoherent  sources,  where  higher  SNR’s  are  required 
in  the  absence  of  noise  bias  compensation. 


B.  Multiple  Signal  Classification  (MUSIC)  Method 

The  MUSIC  method  is  a  second  generation  spectral  estimation  technique  that  primarily 
consists  of  performing  eigenstructure  analysis  on  the  spatial  correlation  matrix  R.  Specifically, 
there  are  two  basic  concepts  behind  the  MUSIC  method  resulting  directly  from  eigenstructure 
analysis.  First,  the  number  of  sources  present  is  determined  from  the  multiplicity  of  the  smallest 
eigenvalue  of  Rg.  Second,  the  directions  of  arrival  are  determined  as  those  whose  direction 
vectors  are  orthogonal  to  the  eigenvectors  corresponding  to  the  smallest  eigenvalue  of  Rj.  Thus, 
the  first  step  of  the  MUSIC  method  is  to  determine  the  complex  eigenvalues  and  corresponding 
eigenvectors  of  R.  The  eigenvectors  are  then  arranged  as  columns  of  an  array  in  accordance  to 
increasing  eigenvalues,  and  the  sum  of  the  outer  products  is  computed  as, 

M 

(25)  Cgy  =  ^  VjVj 

i=l 

where  Vj  is  the  array  of  eigenvectors,  L  is  the  order  of  the  eigenvector  array,  and  q  is  the 
maximum  number  of  sources  that  can  be  detected  as  determined  by  the  multiplicity  of  the  smallest 
eigenvalue  of  R.  Examining  equation  (25)  and  noting  the  arrangement  of  V|,  it  is  apparent  that  the 
summation  is  performed  over  only  those  eigenvectors  that  are  expected  to  correspond  to  signals.- 
Consequently,  the  signals  are  separated  fi’om  the  noise  bias  that  remains  in  the  standard  ME 
method.  Now  the  orthogonality  principles  may  be  applied  to  Cev>  ^^d  the  power  spectrum 
resulting  fi'om  the  MUSIC  method  is  determined  to  be 
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(26) 


SMU=[^WC;vCEvW(u)r‘ 


where  W(u)  is  the  direction  of  look  vector  (equation  (24)).  Schmidt[9]  has  compared  the  MUSIC 
technique  to  the  standard  ME  method  and  the  results  indicate  that  the  MUSIC  method  is  superior. 
This  is  expected  since  the  noise  bias  is  compensated  for  in  the  procedure.  Nevertheless,  the 
MUSIC  approach  is  not  without  limitations.  Consider  the  covariance  or  spatial  correlation  matrix 
described  by  equation  (20).  In  the  absence  of  noise,  R  =  Rg  is  diagonal  when  the  sources  are 
uncorrelated,  non-diagonal  and  non-singular  when  the  sources  are  partially  correlated,  and  non¬ 
diagonal  and  singular  when  some  of  the  sources  are  fully  correlated  or  coherent.  From  basic  linear 
algebra  it  is  apparent  that  this  approach  will  not  perform  well  in  the  presence  of  coherent  sources. 
Consequendy,  matrix  techniques  that  perform  spatial  smoothing  on  R  have  been  developed  to 
remove  the  singularities[10].  These  techniques  are  successful  in  that  when  the  MUSIC  method  is 
applied  to  the  smoothed  spatial  correlation  matrix,  the  performance  in  the  coherent  case  approaches 
that  of  the  incoherent  case.  There  is  however  a  penalty  paid  for  having  to  perform  smoothing. 
Spatial  smoothing  consists  of  diagonal  sub-array  averaging  resulting  in  a  reduction  in  the 
dimensions  of  R,  and  consequently  a  reduction  in  the  number  of  sources  that  may  be  detected. 

C.  Performance  of  Autocorrelation  Based  Techniques 

Up  until  now  the  performance  of  the  ME  and  MUSIC  methods  have  been  discussed  with 
respect  to  a  priori  knowledge  of  the  true  autocorrelation  function.  The  performance  of  all 
autocorrelation  based  techniques  such  as  these,  may  be  characterized  by  recognizing  that  the 
autocorrelation  function  is  not  directly  measured[6].  In  practice,  it  is  estimated  from  the  measmed 
field.  Referring  to  equations  (1)  and  (3),  the  normalized  autocorrelation  function  between  any  two 
elements  spaced  id  apart  is  estimated  by 


N-i 


(27) 


R(ld)  =  x*(n)x(n-^i)  i=0, 1,2,. 

N-i 


n=l 


For  a  two  source  case,  the  sampled  values  of  the  true  autocorrelation  are 
(28)  R(id)  =  expljkiduj]  +  exp[jkidu2]  i=0,l,2,.... 

where  Uj  and  U2  are  the  angular  positions  of  the  point  sources  in  the  far  field.  Using  equations 
(27)  and  (28),  references  [6]  and  [7]  have  thoroughly  evaluated  the  ME  method  in  the  presence  of 
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two  sources.  The  following  discussion  extends  their  results  to  all  autocorrelation  based 
techniques. 

For  the  noise  free  case  of  two  equal  strength,  unit  amplitude,  coherent  point  sources  in  the 
far  field,  the  magnitude  of  the  error  in  the  estimated  autocorrelation  function  is 


(29) 


R(id)-R(id) 


2 

N-i 


sin  kd(N-i) 


sin 


cos(<t)) 


where  ^  is  the  phase  difference  referenced  to  the  center  of  the  array.  Upon  examining  equation 
(29),  some  interesting  conclusions  can  be  made.  Considering  only  the  cosine  term,  the  error  is 
zero  for  (j)  equal  to  odd  multiples  of  7t/2,  and  maximum  for  ())  equal  to  even  multiples  of  7c/2. 
Considering  only  the  ratio  of  the  sine  terms,  as  the  angular  separation  of  the  sources  approaches 
zero,  the  error  approaches  2cos((l) ).  Thus,  the  overall  worst  case  error  will  be  e(id)=2,  which  is 
on  the  order  of  R(0).  Note  that  no  improvement  is  achieved  through  temporal  averaging  in  the 
noiseless  coherent  case;  the  performance  is  identical  regardless  of  the  number  of  snapshots  used. 

If  noise  is  included  in  equation  (29)  the  error  will  become  worse.  Temporal  averaging  can 
help  reduce  the  effects  of  noise,  but  cannot  eliminate  it  (this  will  be  discussed  furtho-  when 
incoherent  sources  are  considered).  Consequently,  if  the  true  autocorrelation  is  unknown,  it 
appears  that  autocorrelation  based  methods  do  not  afford  overall  better  performance  than  that  of  the 
classical  Fourier  method  in  the  case  of  coherent  sources . 

For  the  incoherent  case,  where  the  phase  differences  between  the  sources  varies  with  time 
over  (0,27t),  temporal  averaging  may  be  performed  to  reduce  the  error  described  by  equation  (29) 
to  zero.  Thus,  if  the  duration  of  the  averaging  interval  is  long  enough,  it  is  possible  to  obtain  the 
true  values  of  the  autocorrelation  function.  In  this  case  the  performance  of  the  autocorrelation 
based  methods  are  not  all  identical,  and  for  the  ME  and  MUSIC  approaches,  have  been  previously 
discussed  in  sections  A.  and  B.,  respectively. 
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V.  Conclusions 


The  basic  concepts  behind  spectral  estimation  have  been  presented.  The  classical  Fourier 
approach  exhibits  some  implementation  advantages  that  make  it  a  practical  choice  for  many  low 
performance  applications.  Alternatively,  the  superresolution  techniques  are  computationally 
intensive,  but  can  offer  better  performance  provided  the  limitations  are  understood.  For  the  case  of 
incoherent  sources,  the  ME  method  with  noise  bias  correction  and  the  MUSIC  method  seem  to  be 
comparable,  and  both  offer  better  performance  than  the  classical  Fourier  approach.  For  the  case  of 
coherent  sources,  direct  comparison  is  a  bit  more  involved.  Provided  a  priori  knowledge  of  the 
tme  autocorrelation  function  is  available,  the  ME  method  is  comparable  to  the  classical  Fourier 
approach,  and  both  are  inferior  to  the  MUSIC  method  when  it  is  applied  to  a  smoothed  covariance 
matrix .  If  the  true  autocorrelation  function  is  unknown,  autocorrelation  based  methods  such  as 
ME'and  MUSIC,  do  not  appear  to  afford  overall  better  performance  than  that  of  the  classical 
Fourier  method  for  the  case  of  coherent  sources. 
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